%This script implements the test accounting for estimation uncertainty that
%is reported in Tables 1 and C.1
%Corresponding author: Rodrigo Adao
%Date: 09/11/2024
%Input: model_data_fgkk.mat, est_shares, output_estimation_step.mat
%Output: est_parest.mat, tab_st0

%% Preliminaries
clear;
close all; 

%Define vector of large varieties in our sample
trade_target = 0.9;

%set lacal paths
function_path = "..\Functions\";
addpath(function_path,'-frozen');
set_local_paths;

%cluster: 1 - no cluster, 2 - sector-country 
test_cluster_grid  = 1;

Nspec = length(test_cluster_grid);
results_dW = zeros(17,3,Nspec);
results_dM = zeros(17,3,Nspec);
results_dT = zeros(17,3,Nspec);
results_dX = zeros(17,3,Nspec);

%% Estimation

for j=1:Nspec

tic
    
%% Import and adjust data

display('------------- running new j -----------------')
    j
    
%% Import data

display('Implement preliminary steps')

%estimation parameters    
test_cluster = test_cluster_grid(j);

%Run preliminary steps for estimation
estimation_preliminary_steps;

%% Organzie data from estimation steps

display('Organize data from estimation step')

 %Crosswalk from varieties used in testing to those used in estimation 
    load(save_estimation_path + "output_estimation_step.mat")
    
    %import shifters
    loc_estvariety_Mshifts_mat = zeros(N,G);
    for c = 1:length(all_m_variety_tab)
        loc_estvariety_Mshifts_mat(all_m_variety_tab(c,3),all_m_variety_tab(c,4)) = c;
    end
    
    loc_estvariety_Mshifts = reshape( (  loc_estvariety_Mshifts_mat  )', N*G,1);
    loc_estvariety_Mshifts = loc_estvariety_Mshifts(ind_Mshifts_vec);
    
    %export shifters
    loc_estvariety_Xshifts_mat = zeros(N,G);
    for c = 1:length(all_x_variety_tab)
        loc_estvariety_Xshifts_mat(all_x_variety_tab(c,3),all_x_variety_tab(c,4)) = c;
    end
    
    loc_estvariety_Xshifts = reshape( (  loc_estvariety_Xshifts_mat  )', N*G,1);
    loc_estvariety_Xshifts = loc_estvariety_Xshifts(ind_Xshifts_vec);    

%Organize estimation output for testing
    NTm = 28;
    shiftersM_est = zeros(NMs, NTm);
    shiftersMdiff_est = zeros(NMs, NTm);
    LambdaM_sigma_est = zeros(NMs, NTm);
    LambdaM_omegas_est = zeros(NMs, NTm);
    LambdaM_eta_est = zeros(NMs, NTm);
    LambdaM_kappa_est = zeros(NMs, NTm);

    for c=1:NMs
        if loc_estvariety_Mshifts(c) > 0
           shiftersM_est(c,:)      = delta_tau_est(loc_estvariety_Mshifts(c),:);
           shiftersMdiff_est(c,:)   = diff_tau_est(loc_estvariety_Mshifts(c),:);
           LambdaM_sigma_est(c,:)  = Lambda_sigma(loc_estvariety_Mshifts(c),:);
           LambdaM_omegas_est(c,:) = Lambda_omegas(loc_estvariety_Mshifts(c),:);
           LambdaM_eta_est(c,:)    = Lambda_eta(loc_estvariety_Mshifts(c),:);
           LambdaM_kappa_est(c,:)  = Lambda_kappa(loc_estvariety_Mshifts(c),:);
        end
    end
    Lambda_sigma_est  = [LambdaM_sigma_est ; zeros(NXs,NTm)];
    Lambda_omegas_est = [LambdaM_omegas_est; zeros(NXs,NTm)];
    Lambda_eta_est    = [LambdaM_eta_est   ; zeros(NXs,NTm)];
    Lambda_kappa_est  = [LambdaM_kappa_est ; zeros(NXs,NTm)];

    shiftersX_est = zeros(NXs, NTm);
    shiftersXdiff_est = zeros(NXs, NTm);
    LambdaX_sigmas_est = zeros(NXs, NTm);
    for c=1:NXs
        if loc_estvariety_Xshifts(c) > 0
           shiftersX_est(c,:) = delta_tau_star_est(loc_estvariety_Xshifts(c),:);
           shiftersXdiff_est(c,:) = diff_tau_star_est(loc_estvariety_Xshifts(c),:);
           LambdaX_sigmas_est(c,:) = Lambda_sigmas(loc_estvariety_Xshifts(c),:);
        end
    end
    Lambda_sigmas_est = [zeros(NMs,NTm); LambdaX_sigmas_est];
    
    shifters_est = [shiftersM_est; shiftersX_est];
    shifters_diff_est = [shiftersMdiff_est; shiftersXdiff_est];
    
    Lambda_est = [];
    Lambda_est(:,:,1) = Lambda_sigma_est;
    Lambda_est(:,:,2) = Lambda_omegas_est; 
    Lambda_est(:,:,3) = Lambda_sigmas_est;
    Lambda_est(:,:,4) = Lambda_eta_est;
    Lambda_est(:,:,5) = Lambda_kappa_est;
    
    clearvars Lambda* -except Lambda_est

    %Selecty subset of varieties used in testing with estimation shifters that have at least one tariff change
    ind_shifters_est = ( sum(abs(shifters_diff_est),2) > 0 );
    shifters_est = shifters_est(ind_shifters_est,:);
    shifters_diff_est = shifters_diff_est(ind_shifters_est,:);
    Lambda_est = Lambda_est(ind_shifters_est,:,:);
    
    check_shifters_est = corr([sum(shifters_est,2), sum(shifters_diff_est,2), shiftersIV(ind_shifters_est)]);
    
    %Create shifter IV vector broken down across months for joint estimation
    %wlog we can consider only shifters used in estimation with non-zero changes
    Nm_shifters_est = sum(shifters_diff_est > 0, 2);
    adj_shift_size = full( min( max( shifters(ind_shifters_est)./sum(shifters_diff_est,2) , 0) , 6) ); 
    shifters_est = ( shifters_diff_est.*adj_shift_size - avg_tariff*( (shifters_diff_est > 0)./Nm_shifters_est  ) )/std_tariff;

%Cluster of shifters used for testing and estimation
    if test_cluster == 1
        cluster_est_vec = []; 
    else
        cluster_est_vecM_aux = zeros(length(all_m_variety_tab),1);
        cluster_est_vecM = zeros(NMs,1);
        for c=1:NMs
            if loc_estvariety_Mshifts(c) > 0
               cluster_est_vecM_aux(loc_estvariety_Mshifts(c)) = cluster_vecM(c);
            end
        end
        for c=1:NMs
            if loc_estvariety_Mshifts(c) > 0
               cluster_est_vecM(c) = cluster_est_vecM_aux(loc_estvariety_Mshifts(c));
            end
        end

        cluster_est_vecX_aux = zeros(length(all_x_variety_tab),1);
        cluster_est_vecX = zeros(NXs,1);
        for c=1:NXs
            if loc_estvariety_Xshifts(c) > 0
               cluster_est_vecX_aux(loc_estvariety_Xshifts(c)) = cluster_vecX(c);
            end
        end
        for c=1:NXs
            if loc_estvariety_Xshifts(c) > 0
               cluster_est_vecX(c) = cluster_est_vecX_aux(loc_estvariety_Xshifts(c));
            end
        end

        cluster_est_vec =[cluster_est_vecM; cluster_est_vecX ];
        cluster_est_vec = cluster_est_vec(ind_shifters_est,:);
    end
    clearvars loc_* cluster* -except cluster_est_vec cluster_vec

    %% Compute jacobian of predictions wrt parameters
    display('Compute jacobian of predictions wrt parameters')

    P=5;
    eps=0.01;
    sigma_h      = sigma     *ones(1,P) ;
    omega_star_h = omega_star*ones(1,P);
    sigma_star_h = sigma_star*ones(1,P) ;
    eta_h        = eta       *ones(1,P) ;
    kappa_h      = kappa     *ones(1,P) ;

    sigma_h(1)      = sigma_h(1)      + eps;
    omega_star_h(2) = omega_star_h(2) + eps;
    sigma_star_h(3) = sigma_star_h(3) + eps;
    eta_h(4)        = eta_h(4)        + eps;
    kappa_h(5)      = kappa_h(5)      + eps;

    shareMOD_dW_dtheta = zeros([size(shareMOD_dW), P]);
    for p=1:P
        
        [delta_ig_h, a_star_ig_h, z_star_ig_h, a_ig_h, aM_g_h, AM_s_h, AD_s_h, betaT_h, beta_s_h, alphaL_s_h, alphaI_s_h, alpha_ks_h, barL_rs_h, Z_rs_h, D_h, E_s_h, F_h] = ... 
        invert_param_full(x_ig_0, m_ig_0, tau_star_ig_0, tau_ig_0, trade_con_share, tot_labor_comp, tot_intermediate, sales_s, IO_sales, Lsr, ...
        Dsg, DX, DM, omega_star_h(p), sigma_h(p), eta, kappa, sigma_star_h(p), gammaM, gammaX,  tol_parm);

        
        [p_s_IV, PM_s_IV, ...
                 rho_qm_ig_tau_ig, rho_qm_ig_taustar_ig, rho_qm_ig_a_ig, rho_qm_ig_zstar_ig, rho_qm_ig_astar_ig, rho_qm_ig_delta_ig, ...
                 rho_pm_ig_tau_ig, rho_pm_ig_taustar_ig, rho_pm_ig_a_ig, rho_pm_ig_zstar_ig, rho_pm_ig_astar_ig, rho_pm_ig_delta_ig, ...
                 rho_px_ig_tau_ig, rho_px_ig_taustar_ig, rho_px_ig_a_ig, rho_px_ig_zstar_ig, rho_px_ig_astar_ig, rho_px_ig_delta_ig, ...
                 rho_pstar_ig_tau_ig, rho_pstar_ig_taustar_ig, rho_pstar_ig_a_ig, rho_pstar_ig_zstar_ig, rho_pstar_ig_astar_ig, rho_pstar_ig_delta_ig, ...
                 rho_rm_ig_tau_ig, rho_rm_ig_taustar_ig, rho_rm_ig_a_ig, rho_rm_ig_zstar_ig, rho_rm_ig_astar_ig, rho_rm_ig_delta_ig] ...
                 = compute_equilibrium_foa_full( delta_ig_h, a_star_ig_h, z_star_ig_h, a_ig_h, aM_g_h, AM_s_h, AD_s_h, betaT_h, beta_s_h, alphaL_s_h, alphaI_s_h, alpha_ks_h, barL_rs_h, Z_rs_h, D_h, ... 
                 tau_ig_0, tau_star_ig_0, m_ig_0, E_s_0, p_s_0, Dsg, DX, DM, omega_star_h(p), sigma_h(p), eta_h(p), kappa_h(p), sigma_star_h(p), nu, epsilon, gammaQ, gammaX, gammaM, rhoX, rhoM, tol_conv, adj_inner_g0, adj_outter_g0, ... 
                 ind_Mshifts_vec, ind_Xshifts_vec, ind_Msample_ig, ind_Xsample_ig );  


        share_pm    = [rho_pm_ig_tau_ig    , rho_pm_ig_taustar_ig   ];
        share_px    = [rho_px_ig_tau_ig    , rho_px_ig_taustar_ig   ];
        share_rm    = [rho_rm_ig_tau_ig    , rho_rm_ig_taustar_ig   ];

        shareMOD_dW_dtheta(:,:,p) = (sparse([share_pm; share_rm; share_px]) - shareMOD_dW)/eps;
    end
    

clear    *_h rho_* share_pm share_px share_rm delta_ig_0 a_star_ig_0 z_star_ig_0 a_ig_0 aM_g_0 AM_s_0 AD_s_0 betaT_0 beta_s_0 alphaL_s_0 alphaI_s_0 alpha_ks_0 barL_rs_0 Z_rs_0 ...
             m_ig_0 x_ig_0 tau_ig_0 tau_star_ig_0 Lsr Dsg m_val_sample m_val_adj adj_m_s adj_m_gs x_val_sample x_val_adj adj_x_s adj_x_gs
         
%% IV for welfare

 %Naive IV matrices
    var_adj=sum(ww_dW.^2)/Nobs;
    
    %No GE
    share_NC_nGE_dW  = shareIV_tau;    
    share_wNC_nGE_dW  = ww_dWmat_naive*share_NC_nGE_dW ;
    coef = MC_term1*share_wNC_nGE_dW;
    share_wMC_nGE_dW  = share_wNC_nGE_dW- C*coef;
    share_wMC_nGE_dW  = share_wMC_nGE_dW/var_adj;
    
%Auxiliary matrices for control and welfare adjustments
%load auxiliary adjustment vectors
load(save_share_path + "estimation_shares_t" + trade_target + ".mat", 'adjGE_wMC')

    adjGE_wMC_mat= spdiags(adjGE_wMC, 0 , Nobs, Nobs );
    aux = adjGE_wMC_mat*shareMOD_dW;
    share_wMC_dW_dW = aux - C*(MC_term1*aux);
    
    clearvars shareIV_tau aux adjGE_MC_mat coef ww_dWmat ww_dWmat_naive share_NC_nGE_dW share_wNC_nGE_dW adjGE_wMC_mat

%% Implement test

%Predictions 
    dx = shareMOD_dW*shifters;

    coef = MC_term1*dx;
    dx_MC_dW = dx - C*coef;
    dx_MC_dM = dx_MC_dW(indM==1);
    dx_MC_dT = dx_MC_dW(indT==1);
    dx_MC_dX = dx_MC_dW(indX==1);

    coef = MC_term1*dx(sample_naive);
    dxn_MC_dW = dx(sample_naive) - C*coef;
    dxn_MC_dM = dxn_MC_dW(indMn==1);
    dxn_MC_dT = dxn_MC_dW(indTn==1);
    dxn_MC_dX = dxn_MC_dW(indXn==1);    
    
%Outcomes
    dy = dWout_sample;

    coef = MC_term1*dy;
    dy_MC_dW = dy - C*coef;
    dy_MC_dM = dy_MC_dW(indM==1);
    dy_MC_dT = dy_MC_dW(indT==1);
    dy_MC_dX = dy_MC_dW(indX==1);

    coef = MC_term1*dy(sample_naive);
    dyn_MC_dW = dy(sample_naive) - C*coef;
    dyn_MC_dM = dyn_MC_dW(indMn==1);
    dyn_MC_dT = dyn_MC_dW(indTn==1);
    dyn_MC_dX = dyn_MC_dW(indXn==1);    
    
%Difference between outcome and predictions    
    dyt = dy - dx;

    coef = MC_term1*dyt;
    dyt_MC_dW = dyt - C*coef;
    dyt_MC_dM = dyt_MC_dW(indM==1);
    dyt_MC_dT = dyt_MC_dW(indT==1);
    dyt_MC_dX = dyt_MC_dW(indX==1);

    coef = MC_term1*dyt(sample_naive);
    dytn_MC_dW = dyt(sample_naive) - C*coef;
    dytn_MC_dM = dytn_MC_dW(indMn==1);
    dytn_MC_dT = dytn_MC_dW(indTn==1);
    dytn_MC_dX = dytn_MC_dW(indXn==1);    
    
%jacobian of dy - dx wrt parameters
for p=1:P
    dx_dtheta(:,p) = shareMOD_dW_dtheta(:,:,p)*shifters;

    dyt_dtheta(:,p) = - dx_dtheta(:,p);

    dyt_NC_dW_dtheta(:,p) = dyt_dtheta(:,p);
    coef = MC_term1*dyt_dtheta(:,p);
    dyt_MC_dW_dtheta(:,p) = dyt_dtheta(:,p) - C*coef;

    dyt_MC_dM_dtheta(:,p) = dyt_MC_dW_dtheta(indM==1,p);
    dyt_MC_dT_dtheta(:,p) = dyt_MC_dW_dtheta(indT==1,p);
    dyt_MC_dX_dtheta(:,p) = dyt_MC_dW_dtheta(indX==1,p);

    dytn_NC_dW_dtheta(:,p) = dyt_dtheta(sample_naive,p);
    coef = MC_term1*dyt_dtheta(sample_naive,p);
    dytn_MC_dW_dtheta(:,p) = dyt_dtheta(sample_naive,p) - C*coef;

    dytn_MC_dM_dtheta(:,p) = dytn_MC_dW_dtheta(indMn==1,p);
    dytn_MC_dT_dtheta(:,p) = dytn_MC_dW_dtheta(indTn==1,p);
    dytn_MC_dX_dtheta(:,p) = dytn_MC_dW_dtheta(indXn==1,p);      
end    
    
clear C MC_term1 shareMOD_dW shareMOD_dW_dtheta coef
    
%% Welfare stacked specification

disp('running estimation for welfare stacked outcomes')

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dW, dx_MC_dW]);
    cor_MC_dW = correl_pred(2,1);
    
%Test based on naive, no GE IV
    zn_wMC_dW = share_wMC_nGE_dW*shiftersIV;  
    
    [btt_zn_wMC_dW , set_zn_wMC_dW , rjt_zn_wMC_dW , r0t_zn_wMC_dW, pjt_zn_wMC_dW , p0t_zn_wMC_dW, int_zn_wMC_dW , i0t_zn_wMC_dW ] = implement_jointtest_app_cluster(dytn_MC_dW, zn_wMC_dW , share_wMC_nGE_dW , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dytn_MC_dW_dtheta, ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );          
    [bty_zn_wMC_dW , sey_zn_wMC_dW , rjy_zn_wMC_dW , r0y_zn_wMC_dW, pjy_zn_wMC_dW , p0y_zn_wMC_dW ] = implement_test_cluster(dyn_MC_dW , zn_wMC_dW , share_wMC_nGE_dW , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dW , sex_zn_wMC_dW , rjx_zn_wMC_dW , r0x_zn_wMC_dW, pjx_zn_wMC_dW , p0x_zn_wMC_dW ] = implement_test_cluster(dxn_MC_dW , zn_wMC_dW , share_wMC_nGE_dW , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    zw_wMC_dW = share_wMC_dW_dW*shiftersIV;  
    
    [btt_zw_wMC_dW , set_zw_wMC_dW , rjt_zw_wMC_dW , r0t_zw_wMC_dW, pjt_zw_wMC_dW , p0t_zw_wMC_dW , int_zw_wMC_dW , i0t_zw_wMC_dW] = implement_jointtest_app_cluster(dyt_MC_dW , zw_wMC_dW , share_wMC_dW_dW , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dyt_MC_dW_dtheta , ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );        
    [bty_zw_wMC_dW , sey_zw_wMC_dW , rjy_zw_wMC_dW , r0y_zw_wMC_dW, pjy_zw_wMC_dW , p0y_zw_wMC_dW ] = implement_test_cluster(dy_MC_dW , zw_wMC_dW , share_wMC_dW_dW , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dW , sex_zw_wMC_dW , rjx_zw_wMC_dW , r0x_zw_wMC_dW, pjx_zw_wMC_dW , p0x_zw_wMC_dW ] = implement_test_cluster(dx_MC_dW , zw_wMC_dW , share_wMC_dW_dW , shiftersIV, critical, 0, cluster_vec);

%Output
    results_zw_dW = full([bty_zw_wMC_dW, btx_zw_wMC_dW, btt_zw_wMC_dW; ...
                          sey_zw_wMC_dW, sex_zw_wMC_dW, set_zw_wMC_dW; ...
                          pjy_zw_wMC_dW, pjx_zw_wMC_dW, pjt_zw_wMC_dW; ...
                          p0y_zw_wMC_dW, p0x_zw_wMC_dW, p0t_zw_wMC_dW; ...
                                      0,             0, int_zw_wMC_dW ; ...
                                      0,             0, i0t_zw_wMC_dW]);

    results_zn_dW = full([bty_zn_wMC_dW, btx_zn_wMC_dW, btt_zn_wMC_dW; ...
                          sey_zn_wMC_dW, sex_zn_wMC_dW, set_zn_wMC_dW; ...
                          pjy_zn_wMC_dW, pjx_zn_wMC_dW, pjt_zn_wMC_dW; ...
                          p0y_zn_wMC_dW, p0x_zn_wMC_dW, p0t_zn_wMC_dW; ...
                                      0,             0, int_zn_wMC_dW ; ...
                                      0,             0, i0t_zn_wMC_dW]);

    results_cor_dW =full([0             , 0             , cor_MC_dW    ]);

    results_obs_dW = Nobs*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dW_j = [results_zw_dW; bar; results_zn_dW; bar; results_cor_dW; bar; results_obs_dW]; 
    tab_outcomes_dW = table(indX, indM, indT, ww_dW, dy_MC_dW, dx_MC_dW, dyt_MC_dW, zn_wMC_dW, zw_wMC_dW) ;

    clear results_zw_dW results_zn_dW results_cor_dW results_obs_dW  bt* se* rj* r0* pj* p0*
    
%% Import specification

disp('running estimation for import outcomes') 

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dM, dx_MC_dM]);
    cor_MC_dM = correl_pred(2,1);
    
%Test based on naive, no GE IV
    share_wMC_nGE_dM  = share_wMC_nGE_dW(indMn==1,:);
    zn_wMC_dM = share_wMC_nGE_dM*shiftersIV;  
    
    [btt_zn_wMC_dM , set_zn_wMC_dM , rjt_zn_wMC_dM , r0t_zn_wMC_dM, pjt_zn_wMC_dM , p0t_zn_wMC_dM, int_zn_wMC_dM , i0t_zn_wMC_dM ] = implement_jointtest_app_cluster(dytn_MC_dM, zn_wMC_dM , share_wMC_nGE_dM , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dytn_MC_dM_dtheta, ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );          
    [bty_zn_wMC_dM , sey_zn_wMC_dM , rjy_zn_wMC_dM , r0y_zn_wMC_dM, pjy_zn_wMC_dM , p0y_zn_wMC_dM ] = implement_test_cluster(dyn_MC_dM , zn_wMC_dM , share_wMC_nGE_dM , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dM , sex_zn_wMC_dM , rjx_zn_wMC_dM , r0x_zn_wMC_dM, pjx_zn_wMC_dM , p0x_zn_wMC_dM ] = implement_test_cluster(dxn_MC_dM , zn_wMC_dM , share_wMC_nGE_dM , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    share_wMC_dW_dM  = share_wMC_dW_dW(indM==1,:);
    zw_wMC_dM = share_wMC_dW_dM*shiftersIV;  
    
    [btt_zw_wMC_dM , set_zw_wMC_dM , rjt_zw_wMC_dM , r0t_zw_wMC_dM, pjt_zw_wMC_dM , p0t_zw_wMC_dM , int_zw_wMC_dM , i0t_zw_wMC_dM] = implement_jointtest_app_cluster(dyt_MC_dM , zw_wMC_dM , share_wMC_dW_dM , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dyt_MC_dM_dtheta , ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );        
    [bty_zw_wMC_dM , sey_zw_wMC_dM , rjy_zw_wMC_dM , r0y_zw_wMC_dM, pjy_zw_wMC_dM , p0y_zw_wMC_dM ] = implement_test_cluster(dy_MC_dM , zw_wMC_dM , share_wMC_dW_dM , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dM , sex_zw_wMC_dM , rjx_zw_wMC_dM , r0x_zw_wMC_dM, pjx_zw_wMC_dM , p0x_zw_wMC_dM ] = implement_test_cluster(dx_MC_dM , zw_wMC_dM , share_wMC_dW_dM , shiftersIV, critical, 0, cluster_vec);

%Output
    results_zw_dM = full([bty_zw_wMC_dM, btx_zw_wMC_dM, btt_zw_wMC_dM; ...
                          sey_zw_wMC_dM, sex_zw_wMC_dM, set_zw_wMC_dM; ...
                          pjy_zw_wMC_dM, pjx_zw_wMC_dM, pjt_zw_wMC_dM; ...
                          p0y_zw_wMC_dM, p0x_zw_wMC_dM, p0t_zw_wMC_dM; ...
                                      0,             0, int_zw_wMC_dM ; ...
                                      0,             0, i0t_zw_wMC_dM]);

    results_zn_dM = full([bty_zn_wMC_dM, btx_zn_wMC_dM, btt_zn_wMC_dM; ...
                          sey_zn_wMC_dM, sex_zn_wMC_dM, set_zn_wMC_dM; ...
                          pjy_zn_wMC_dM, pjx_zn_wMC_dM, pjt_zn_wMC_dM; ...
                          p0y_zn_wMC_dM, p0x_zn_wMC_dM, p0t_zn_wMC_dM; ...
                                      0,             0, int_zn_wMC_dM ; ...
                                      0,             0, i0t_zn_wMC_dM]);

    results_cor_dM =full([0             , 0             , cor_MC_dM    ]);

    results_obs_dM = NM*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dM_j = [results_zw_dM; bar; results_zn_dM; bar; results_cor_dM; bar; results_obs_dM]; 

    clear results_zw_dM results_zg_dM results_zn_dM results_cor_dM results_obs_dM share_wMC_nGE_dM share_wMC_dW_dM bt* se* rj* r0* pj* p0* 
    
 %% Revenue specification

disp('running estimation for revenue outcomes')

%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dT, dx_MC_dT]);
    cor_MC_dT = correl_pred(2,1);
    
%Test based on naive, no GE IV
    share_wMC_nGE_dT  = share_wMC_nGE_dW(indTn==1,:);
    zn_wMC_dT = share_wMC_nGE_dT*shiftersIV;  
    
    [btt_zn_wMC_dT , set_zn_wMC_dT , rjt_zn_wMC_dT , r0t_zn_wMC_dT, pjt_zn_wMC_dT , p0t_zn_wMC_dT, int_zn_wMC_dT , i0t_zn_wMC_dT ] = implement_jointtest_app_cluster(dytn_MC_dT, zn_wMC_dT , share_wMC_nGE_dT , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dytn_MC_dT_dtheta, ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );          
    [bty_zn_wMC_dT , sey_zn_wMC_dT , rjy_zn_wMC_dT , r0y_zn_wMC_dT, pjy_zn_wMC_dT , p0y_zn_wMC_dT ] = implement_test_cluster(dyn_MC_dT , zn_wMC_dT , share_wMC_nGE_dT , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dT , sex_zn_wMC_dT , rjx_zn_wMC_dT , r0x_zn_wMC_dT, pjx_zn_wMC_dT , p0x_zn_wMC_dT ] = implement_test_cluster(dxn_MC_dT , zn_wMC_dT , share_wMC_nGE_dT , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    share_wMC_dW_dT  = share_wMC_dW_dW(indT==1,:);
    zw_wMC_dT = share_wMC_dW_dT*shiftersIV;  
    
    [btt_zw_wMC_dT , set_zw_wMC_dT , rjt_zw_wMC_dT , r0t_zw_wMC_dT, pjt_zw_wMC_dT , p0t_zw_wMC_dT , int_zw_wMC_dT , i0t_zw_wMC_dT] = implement_jointtest_app_cluster(dyt_MC_dT , zw_wMC_dT , share_wMC_dW_dT , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dyt_MC_dT_dtheta , ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );        
    [bty_zw_wMC_dT , sey_zw_wMC_dT , rjy_zw_wMC_dT , r0y_zw_wMC_dT, pjy_zw_wMC_dT , p0y_zw_wMC_dT ] = implement_test_cluster(dy_MC_dT , zw_wMC_dT , share_wMC_dW_dT , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dT , sex_zw_wMC_dT , rjx_zw_wMC_dT , r0x_zw_wMC_dT, pjx_zw_wMC_dT , p0x_zw_wMC_dT ] = implement_test_cluster(dx_MC_dT , zw_wMC_dT , share_wMC_dW_dT , shiftersIV, critical, 0, cluster_vec);

%Output
    results_zw_dT = full([bty_zw_wMC_dT, btx_zw_wMC_dT, btt_zw_wMC_dT; ...
                          sey_zw_wMC_dT, sex_zw_wMC_dT, set_zw_wMC_dT; ...
                          pjy_zw_wMC_dT, pjx_zw_wMC_dT, pjt_zw_wMC_dT; ...
                          p0y_zw_wMC_dT, p0x_zw_wMC_dT, p0t_zw_wMC_dT; ...
                                      0,             0, int_zw_wMC_dT ; ...
                                      0,             0, i0t_zw_wMC_dT]);

    results_zn_dT = full([bty_zn_wMC_dT, btx_zn_wMC_dT, btt_zn_wMC_dT; ...
                          sey_zn_wMC_dT, sex_zn_wMC_dT, set_zn_wMC_dT; ...
                          pjy_zn_wMC_dT, pjx_zn_wMC_dT, pjt_zn_wMC_dT; ...
                          p0y_zn_wMC_dT, p0x_zn_wMC_dT, p0t_zn_wMC_dT; ...
                                      0,             0, int_zn_wMC_dT ; ...
                                      0,             0, i0t_zn_wMC_dT]);

    results_cor_dT =full([0             , 0             , cor_MC_dT    ]);

    results_obs_dT = NM*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dT_j = [results_zw_dT; bar; results_zn_dT; bar; results_cor_dT; bar; results_obs_dT]; 

    clear results_zw_dT results_zg_dT results_zn_dT results_cor_dT results_obs_dT share_wMC_nGE_dT share_wMC_dW_dT bt* se* rj* r0* pj* p0* 
    
  %% Export specification

disp('running estimation for export outcomes')


%Correlation correlation-based tests
    correl_pred = corr([dy_MC_dX, dx_MC_dX]);
    cor_MC_dX = correl_pred(2,1);
    
%Test based on naive, no GE IV
    share_wMC_nGE_dX  = share_wMC_nGE_dW(indXn==1,:);
    zn_wMC_dX = share_wMC_nGE_dX*shiftersIV;  
    
    [btt_zn_wMC_dX , set_zn_wMC_dX , rjt_zn_wMC_dX , r0t_zn_wMC_dX, pjt_zn_wMC_dX , p0t_zn_wMC_dX, int_zn_wMC_dX , i0t_zn_wMC_dX ] = implement_jointtest_app_cluster(dytn_MC_dX, zn_wMC_dX , share_wMC_nGE_dX , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dytn_MC_dX_dtheta, ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );          
    [bty_zn_wMC_dX , sey_zn_wMC_dX , rjy_zn_wMC_dX , r0y_zn_wMC_dX, pjy_zn_wMC_dX , p0y_zn_wMC_dX ] = implement_test_cluster(dyn_MC_dX , zn_wMC_dX , share_wMC_nGE_dX , shiftersIV, critical, 0, cluster_vec);
    [btx_zn_wMC_dX , sex_zn_wMC_dX , rjx_zn_wMC_dX , r0x_zn_wMC_dX, pjx_zn_wMC_dX , p0x_zn_wMC_dX ] = implement_test_cluster(dxn_MC_dX , zn_wMC_dX , share_wMC_nGE_dX , shiftersIV, critical, 0, cluster_vec);

%Test based on welfare GE IV
    share_wMC_dW_dX  = share_wMC_dW_dW(indX==1,:);
    zw_wMC_dX = share_wMC_dW_dX*shiftersIV;  
    
    [btt_zw_wMC_dX , set_zw_wMC_dX , rjt_zw_wMC_dX , r0t_zw_wMC_dX, pjt_zw_wMC_dX , p0t_zw_wMC_dX , int_zw_wMC_dX , i0t_zw_wMC_dX] = implement_jointtest_app_cluster(dyt_MC_dX , zw_wMC_dX , share_wMC_dW_dX , shiftersIV, shifters_est, Lambda_est, Vgamma, Jgamma, dyt_MC_dX_dtheta , ind_shifters_est, critical , 0, cluster_vec, cluster_est_vec );        
    [bty_zw_wMC_dX , sey_zw_wMC_dX , rjy_zw_wMC_dX , r0y_zw_wMC_dX, pjy_zw_wMC_dX , p0y_zw_wMC_dX ] = implement_test_cluster(dy_MC_dX , zw_wMC_dX , share_wMC_dW_dX , shiftersIV, critical, 0, cluster_vec);
    [btx_zw_wMC_dX , sex_zw_wMC_dX , rjx_zw_wMC_dX , r0x_zw_wMC_dX, pjx_zw_wMC_dX , p0x_zw_wMC_dX ] = implement_test_cluster(dx_MC_dX , zw_wMC_dX , share_wMC_dW_dX , shiftersIV, critical, 0, cluster_vec);

%Output
    results_zw_dX = full([bty_zw_wMC_dX, btx_zw_wMC_dX, btt_zw_wMC_dX; ...
                          sey_zw_wMC_dX, sex_zw_wMC_dX, set_zw_wMC_dX; ...
                          pjy_zw_wMC_dX, pjx_zw_wMC_dX, pjt_zw_wMC_dX; ...
                          p0y_zw_wMC_dX, p0x_zw_wMC_dX, p0t_zw_wMC_dX; ...
                                      0,             0, int_zw_wMC_dX ; ...
                                      0,             0, i0t_zw_wMC_dX]);

    results_zn_dX = full([bty_zn_wMC_dX, btx_zn_wMC_dX, btt_zn_wMC_dX; ...
                          sey_zn_wMC_dX, sex_zn_wMC_dX, set_zn_wMC_dX; ...
                          pjy_zn_wMC_dX, pjx_zn_wMC_dX, pjt_zn_wMC_dX; ...
                          p0y_zn_wMC_dX, p0x_zn_wMC_dX, p0t_zn_wMC_dX; ...
                                      0,             0, int_zn_wMC_dX ; ...
                                      0,             0, i0t_zn_wMC_dX]);

    results_cor_dX =full([0             , 0             , cor_MC_dX    ]);

    results_obs_dX = NX*ones(1,3);                     
    
    bar = zeros(1,3);
    results_dX_j = [results_zw_dX; bar; results_zn_dX; bar; results_cor_dX; bar; results_obs_dX]; 

    clear results_zw_dX results_zg_dX results_zn_dX results_cor_dX results_obs_dX share_wMC_nGE_dX share_wMC_dW_dX bt* se* rj* r0* pj* p0* 
    
%% save output

    results_dW(:,:, j) = results_dW_j;
    results_dM(:,:, j) = results_dM_j;
    results_dT(:,:, j) = results_dT_j;
    results_dX(:,:, j) = results_dX_j;

    clear share_wMC_nGE_dW share_wMC_dW_dW zw* zn* *_j *_dtheta

    if j == 1
        tab_st0 = tab_outcomes_dW; 
    end      

    toc
end

save(save_estimation_path + "est_Tab1b_TabC1.mat", 'results_dW', 'results_dM', 'results_dT', 'results_dX', 'test_cluster_grid', 'tab_st0', 'trade_target');
writetable(tab_st0,save_estimation_path + "tab_st0" )

%% Visualizing results

load(save_estimation_path + "est_Tab1b_TabC1.mat")

j=1;
Tab1 = results_dW(:,:, j) ;

Tab1_panelB = [Tab1(15,3), Tab1(8,3) , Tab1(1,3); 
                0        , Tab1(9,3) , Tab1(2,3);
                0        , Tab1(11,3), Tab1(4,3)];

Tab_C1 =    [   results_dW(1:4  ,:, j)  ;
              results_dW(end  ,:, j) ;
              zeros(1,3);
              results_dX(1:4  ,:, j)  ;
              results_dX(end  ,:, j)  ;
              zeros(1,3);
              -results_dM(1  ,:, j)  ;
              results_dM(2:4  ,:, j)  ;
              results_dM(end  ,:, j)  ;
              zeros(1,3);
              results_dT(1:4  ,:, j)  ;
              results_dT(end  ,:, j)  ];
    

